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We study the effects of gravitational lensing on the estimation of non-Gaussianity from the bis- 
pectrum of the cosmic microwave background (CMB) temperature anisotropies. We find that the 
effect of lensing on the bispectrum may qualitatively be described as a smoothing of the acous- 
tic features analogous to the temperature power spectrum. In contrast to previous results, for a 
Planck-like experiment which is cosmic-variance limited to ^ ma x = 2000, we find that lensing causes 
no significant degradation of our ability to constrain the non-Gaussianity amplitude f^L for both 
local and equilateral configurations, provided that the biases due to the cross correlation between 
the lensing potential and the integrated-Sachs- Wolfe (ISW) contribution to the CMB temperature 
are adequately understood. With numerical simulations, we also verify that low-order Taylor ap- 
proximations to the lensed bispectrum and ISW-lensing biases are accurate. 



I. INTRODUCTION 

It has become clear that primordial non-Gaussianity 
is a powerful tool to constrain different models of in- 
flation and shed light on the physics of the early Uni- 
verse. A large range of early-Universe models are com- 
patible with current measurements of the CMB power 
spectrum, provided that they can produce small (nearly) 
scale-invariant primordial curvature perturbations in an 
otherwise flat universe. Distinguishing amongst these 
models will require not only additional measurements, 
but also characterizations of the data beyond the power 
spectrum. The first such statistic which is available is the 
bispectrum or three-point correlation function in Fourier 
space. As current observations already constrain the non- 
Gaussianity of the CMB to be weak, it can be shown that 
the bispectrum is also an optimal statistic to study [T], 
and so it has justifiably become the subject of much work. 

The primordial bispectrum B(ki,k 2 ,k 3 ) is usually 
characterized by an overall amplitude, given by the 
dimensionless parameter Jnl, and a shape specifying 
which configurations of wavevectors contain the highest 
contributions to the non-Gaussian signal. Translational 
invariance imposes the constraint k\ + k 2 + ks = 0, and 
rotational and parity invariance forces the bispectrum to 
be a function of the lengths of the three wavevectors only. 
Thus, bispectrum shapes are often idealized as those of 
triangles. The two most common choices are the local 
shape (hereafter often denoted as loc), where the signal 
is maximum on squeezed configurations (fci <C k 2 ,ks); 
and the equilateral shape (hereafter often denoted as 
eq), in which the bispectrum peaks mostly on equilat- 
eral triangles {k\ ~ k 2 ~ ^3). Many scenarios for the 
generation of the primordial curvature perturbation fall 
more or less into one of these classes. The local shape 
is generally produced by models in which the perturba- 
tions are generated outside the horizon, curvaton models 
[2H4] . In single- field inflation, the local shape cannot be 
generated at a detectable (i.e. > 0(1)) level; there is 



a theorem which states that the single-field bispectrum 
in squeezed triangles is proportional to the tilt (1 — n s ) 
of the power spectrum [5] , and current observations con- 
strain the power spectrum to be nearly scale-invariant 
6J . Equilateral shapes are a signature of nonstandard ki- 
netic terms in the inflaton Lagrangian, as for example in 
DBI [7 and ghost inflation [5] . For a complete discussion 
on shape classification of primordial bispectra and their 
correlations see [9]. This standard classification scheme 
provides a very useful interface between observation and 
theory. It allows analysts to focus on constraining the 
amplitudes of only the fundamental bispectra, and it al- 
lows theorists to check rapidly whether their models are 
consistent with current observational constraints. 

The current best 2a observational limits on Jnl pa- 
rameters from the WMAP 5-year data are [BJ [THI ITTj : 
-4 < ffil < 80 and -125 < /^ q L < 435. Combin- 
ing WMAP and SDSS data pi] yields -K /#£ < 63. 
Thus the current data do not support a detection of non- 
Gaussianity, although the evidence for f£ c L is close to 
2cr. These results will soon improve dramatically: fore- 
casted uncertainties from the Planck satellite are roughly 
~ 5 [13] and cr(/^ q L ) « 60 0. Such an im- 
provement on the present error bars should allow us to 
tighten significantly our present constraints on inflation- 
ary scenarios. A detection of primordial f l £ c L > 1, for 
example, would rule out standard single-field slow-roll 
inflation [T51 HH] — a potentially sea-changing result. 

Given the deep implications that a detection of pri- 
mordial non-Gaussianity would have, it is crucial that all 
possible sources of contamination for the non-Gaussian 
measurement are well under control. In other words, we 
have to make sure that if a signal is extracted from CMB 
data using estimators of non-Gaussianity, it is of primor- 
dial origin and not produced by some spurious secondary 
or instrumental effect. 

Many different sources could in principle bias a pri- 
mordial non-Gaussianity measurement. In the analyses 
of WMAP data performed so far, particular attention 
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has been devoted to astrophysical contaminants such as 
residual foreground contamination and unresolved point 
sources [5J [TO] E] as potential spurious signals. Another 
possible source of contamination is the non-Gaussianity 
induced by second-order anisotropics. Beyond linear or- 
der in perturbation theory, it is no longer true that Gaus- 
sian initial conditions imply Gaussianity of the CMB 
temperature field. It is therefore important to study 
secondary anisotropics that produce non-Gaussianities of 
similar amplitude and shape as the primordial ones in the 
CMB. In order to study this aspect in a fully consistent 
way, a complete numerical implementation of the second- 
order Einstein-Boltzmann evolution equations [T%rt2"T] is 
necessary. Only a partial implementation is available at 
present [Hj. Meanwhile, in the absence of a full numeri- 
cal solution, a number of papers on the subject [HI |2"3T - 
[29] have focused on specific well-known secondaries such 
as gravitational lensing, the Sunyaev-Zel'dovich effect, 
and perturbed recombination. Their effects have been 
found small for £ < 500, and so do not form a significant 
source of contamination for WMAP. For higher resolution 
experiments such as Planck, however, they are expected 
to dominate over e.g. residual point sources, and must 
be treated with care [26] . 

In this paper we focus on the secondary non- 
Gaussianity induced by gravitational lensing of the CMB. 
Note that while lensing itself does not generate a three- 
point function, if the lensing effects are correlated to 
the unlensed CMB then a bispectrum may be generated. 
Such a correlation arises at low-£ from the integrated 
Sachs- Wolfe (ISW) effect or at high-^ from the nonlinear 
ISW (Rees-Sciama) and Sunyaev-Zel'dovich effects. The 
ISW-lcnsing bispectrum is a direct source of bias when 
estimating the /nl parameters, and can "fake" the pri- 
mordial signal if not accounted for in the analysis. This is 
particularly true for the local shape, as the ISW-lensing 
correlation sources squeezed bispectrum modes. Lensing 
is expected to provide the largest source of secondary bias 
for f y £ c L estimation [2"rj] . 

Apart from this direct ISW-lensing bias, the shape of 
the observed bispectrum is also modified by lensing. This 
shape change could modify the effective normalization for 
a /jvl estimator or even confuse the different primordial 
shapes with each other. In [35], for example, it was found 
that lensing generates a large change in the shape of the 
observed three-point function, degrading the experimen- 
tal sensitivity to /jvl (as will be discussed later, we do 
not reproduce this result). 

Both the ISW-lensing bias and shape change due to 
lensing can be approximated analytically. Gravitational 
lensing of the CMB is treated as a deflection of the lines 
of sight between the observer and recombination, with 
preserved surface brightness (for a recent review see [3D]). 
The lensed CMB T(n) is related to the unlensed CMB 
T(n) by 



where <f> is the lensing potential. 1 For analytical calcula- 
tions involving T(n), Eq. is usually Taylor expanded 
to second order in <j>, and ensemble-averaged results are 
taken to O(Cf^). The accuracy of this approximation 
has been studied thoroughly in the context of the lensed 
power spectrum |32j . where it results in errors of order 
10% of the lensing effect at intermediate multipoles. A 
primary purpose of this paper is to verify with simula- 
tions of the exact lensing displacements that a low-order 
Taylor approximation is similarly accurate for the bis- 
pectrum. 

The remainder of this paper is organized as follows. 
In f|TT] we detail our simulation and analysis steps. We 
calculate the bias due to ISW-lensing in §111 A and we 
investigate the effect of lensing on the shape and normal- 
ization of the primordial bispectrum in §111 B[ In §111 C| 
we study the increased statistical error in /jvl parame- 
ters due to non-Gaussian statistics of the lensed CMB. 
We summarize and draw our conclusions in §IV[ The ap- 
pendix provides more details of our simulation method- 
ology. 



II. PRELIMINARIES AND NOTATION 



The angular bispectrum Bg^^ is defined by 



(ae imi a,£ 2m2 ae 3m3 } - B £l i 2 i 3 ( ^ ^ ^ ) . (2) 



Here, the a £m are the spherical-multipole coefficients of 
the observed CMB and the ensemble average is taken over 
realizations of the primordial perturbations. This is the 
most general form of the three-point function which is 
rotationally invariant. Under the additional assumption 
of parity invariance (so that Bi x i 2 i 3 = if l\ + 4 + 4 is 
odd and so it is invariant under all permutations) we can 
define the reduced bispectrum bf,^i 3 by 
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The reduced bispectra for the local and equilateral shapes 
can be computed efficiently as integrals involving the 
CMB transfer functions and the primordial power spec- 
trum; see, for example, [53] for details. The CMB bis- 
pectra are characterized by an amplitude- f$ L , where X 
denotes either local or equilateral, and a shape. We shall 
generally denote the primordial CMB bispectra with unit 
f NL (i.e. the shape part) by Bf ll2 n 3 . 

In the limit of weak non-Gaussianity, the minimum- 
variance full-sky estimator for f^ L given cosmic-variance 



T(n) = T[n + V0(n)] 



(1) 



1 For a discussion of the spherical displacements that are implied 
by Eq. jl}, see Ref. [31]. 
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limited data to 



is given by 



III. RESULTS 
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where the Fisher-matrix element F(B, B') is defined for 
bispectra B, B' by 



(5) 



Note that F~ 1 (B X , B x ) gives the variance of the error 
in in the Gaussian approximation. The harmonic- 
space form of the estimator in Eq. Q is too slow for prac- 
tical use, but there is a mathematically equivalent, fast 
position-space form for the local and equilateral shapes 



We use non-Gaussian CMB simulations both to verify 
the accuracy of our low-order analytical results and to 
approximate quantities which are too intensive to calcu- 
late directly. Our simulations are composed of "pairs" of 
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This is discussed further in the 



The Gaussian part af m sets the power spectrum of our 
simulations, while the non-Gaussian part af^ sets the 
shape of the bispectrum. We generate the non-Gaussian 
part from af m using the simulation algorithm from [14] . 
This approach is designed for maximum generality and 
allows us to construct non-Gaussian simulations for any 
bispectrum, although for specific bispectra it often con- 
tains freedoms which may be adjusted to modify the vari 
ance of the realized a 
appendix. 

Our simulations are of a flat ACDM cosmology with 
{n b ,fl c ,h,n s ,T,A s } = {0.05,0.23,0.7,0.96,0.08,2.4 x 
10~ 9 }. We simulate the lensing potential 0(n) as a Gaus- 
sian field which is correlated to the unlensed CMB via the 
ISW effect. The auto power spectra Cf T , Cg/ and cross 

spectrum Cj* are computed using C AMB [35] . We sim- 
ulate the deflection operation of Eq. ([I]) using the public 
LENSPIX code (351 [37] , which performs cubic interpola- 
tion on a high-resolution map. LENSPIX produces re- 
sults which are consistent with the "true" lensed power 
spectrum calculated following [33] to 0.1% at i < 2000. 

For our /jvl analysis, we will restrict the sum of Eq. Q 
to £ max = 2000 to mimic the cosmic-variance limit ex- 
pected from the Planck satellite. For lensed and unlensed 
simulations, we use lensed and unlensed power spectra 
respectively in the denominator of the estimator |38j . 



A. ISW-lensing bias 

The most worrying source of contamination for an 
analysis of primordial non-Gaussianity is the nonzero 
bispectrum generated by lensing, since this can directly 
bias the /jvl estimators. If the lens potential 4> and the 
unlensed CMB are statistically independent, then lens- 
ing cannot generate a bispectrum, because there is a 
T — > (— T) symmetry. This argument does not make any 
approximations (such as expanding to a finite order in 
powers of the lens potential, or assuming that the lens po- 
tential is a Gaussian field). Because there is an ISW cross 
correlation, however, lensing can generate a bispectrum. 
We are interested in the bias (/jvxhsw-iensing to the /jvl 
estimators. The ISW-lensing bispectrum can also be used 
as a source of cosmological information [211 [Ml SO] > but 
here we concentrate on the bias to the primordial ampli- 
tudes. 

To lowest order in Cj^ , we can easily predict the bias. 
The ISW-lensing bispectrum is 
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The resultant bias to f x L can be obtained by comput- 
ing the expectation value of the estimator in Eq. Q in 
the presence of the approximate ISW-lensing bispectrum 
[Eq. Q]. A short calculation shows that 



/jvl / 

/ ISW-lonsing 



F( B X , 5 (ISW- lensing) \ 
F{B X ,B X ) ' 



(9) 



This bias is plotted in Fig. [T] for both the local and equi- 
lateral shapes as a function of £ mllx for a cosmic- variance 
limited /jvl estimator. For local configurations, the bias 
is significant, but for equilateral configurations it is al- 
ways at most one order of magnitude below the estima- 
tor variance, as has been observed elsewhere [21 |2"6] , 
This behavior follows since large-scale potential fluctu- 
ations source the ISW effect and also lens the CMB on 
small scales, producing a bispectrum in squeezed trian- 
gles which is correlated with the local shape. 

Equation ^ is only a leading-order approximation, 
when the lensing operation in Eq. ([T]) is expanded in pow- 
ers of </>. We can use the lensed simulations described in 
pi] to test the accuracy of this approximation when com- 
puting biases in /jvl ■ Table [I] compares the lowest-order 
approximate results with those from 100 Monte-Carlo 
simulations for an experiment that is cosmic-variance 
limited to £ max = 2000. (The errors quoted for the sim- 
ulation results are the standard error in the mean /jvl 
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FIG. 1: Biases in f§ L for the local (top) and equilateral (bot- 
tom) shapes if the ISW-lensing cross correlation were to be 
ignored. The analysis is assumed cosmic-variance limited up 
to a maximum multipole i? ma x- The solid/dotted lines are cal- 
culated from Eq. |9| and are shown dotted where the bias is 
negative. Long-dashed lines are the expected Gaussian errors 
on f* L computed from the Fisher matrix. 



Fisher Simulations 


Local 


+9.3 +9.4 ±0.2 


Equilateral 


-2.4 -3.1 + 1.8 



TABLE I: Biases in f$ L from the ISW-lensing cross correla- 
tion. 



from 100 simulations.) The agreement is excellent. As 
a technical point, to reduce the measurement error from 
the finite Monte-Carlo sample, we have subtracted the 
spurious contribution to Jnl from the unlensed CMB in 
each Monte-Carlo realization, i.e. we estimate the ISW- 
lensing bias as follows: 



ISW-lensing 



rX [" lensed] fX [" unlensed] 
JNL l a (m \ — JNL Y a tm J 



(10) 

The second term on the right-hand side has zero mean 
since the unlensed CMB has a vanishing three-point func- 
tion, but including it improves the statistical error on 
}isw-iensing due to the finite Monte-Carlo sample. 
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B. Lensing of the primordial bispectrum 

Gravitational lensing can also change the "shape" of 
the primordial bispectrum. This is a milder form of con- 
tamination than the ISW-lensing signal in §111 A| It can- 
not fake a primordial signal, but it can confuse the dif- 
ferent shapes with each other or degrade the sensitivity 



of the estimator. 

To first order in Cf^, the lensed bispectrum Bi x i % i 3 is 
given by [2"5] 

where B^e^ denotes the unlensed bispectrum, and we 
have defined 
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Here, the symbol in braces is a Wigner-6j symbol. It 
is straightforward to show that the lensed bispectrum 
inherits the symmetry properties of the unlensed bispec- 
trum, as required, since lensing respects rotational and 
parity invariance. 



Evaluating Eq. (11) for both the local and equilateral 
shapes, we find that the difference between the lensed 
bispectrum Bg 1 g 2 i 3 and the unlensed bispectrum Bg 1 g 2 g 3 
is small. This disagrees substantially with the findings 
of |25j . although our analytical approach is the same. In 
Fig. [2] we plot the lensing effects on a slice through the 
local bispectrum. The agreement between our simula- 
tions and analytical results is excellent, with only small 
discrepancies at high-£ due to higher-order terms in Cf^ 
that are neglected in the analytic calculation. The main 
effect of lensing on the bispectrum is a smoothing of its 
acoustic features, analogous to the lensing corrections to 
the TT and EE power spectra. The magnitude of the ef- 
fect is also quantitatively very similar to the power spec- 
trum case, on the order of 10%. 

The agreement between the 0(Cf^) result for the 
lensed bispectrum and our simulations is rather better 
than a casual inspection of Eq. (12 1 might suggest. The 



first term arises from three-point correlations of the form 
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((V^Si)) (V^(ni)) (V'VTCnO) T(n 2 )T(n 3 )> 



= -(a 2 )(V 2 T(n 1 )T(n 2 )T(n 3 )), 



(13) 



i.e. the unlensed CMB at two points correlated with the 
second-order term in the Taylor expansion of the lensed 
CMB at a third point. Here, (ex 2 ) ~ (2.7arcmin) 2 is the 
mean-square of the lensing deflection angle (a = V0). 
Including higher-order terms in the Taylor expansion of 
the lensing action, the first term in Eq. ( 12 1 therefore gen- 
eralizes to [exp(-£ 2 (a 2 )/4) - l]B^ t J 3 for I > 1. This 

is poorly approximated by an expansion to 0(Cf ) for 
I > 1000 and is symptomatic of poor converge of the Tay- 
lor expansion of the lensing action in the map on small 
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FIG. 2: A "slice" Be,i+io,io through the local bispectrum for 
}nl = 1. The simulations are unlensed (magenta line) and 
lensed (cyan line) with Cf^ = 0, while the first-order analyt- 
ical predictions are solid black (unlensed) and dashed black 
(lensed) lines. The Monte-Carlo results use fOOO simulations. 
The fractional effect due to lensing is shown in the bottom 
panel for the simulations (red line) and the first-order ana- 
lytic result (black line). 



scales. However, what matters for the statistics of the 
lensed CMB is the relative displacement of points and 
this is much better approximated by a low-order Taylor 
expansion than the absolute displacements. When calcu- 
lating the lensed A^-point functions in Fourier space, the 
importance of only relative displacements is hidden in a 
near cancellation between (large) terms of the sa me o rder 
in C//. For example, the second term in Eq. (12 1 and 



its cyclic permutations, which arise from correlations of 
the form 



(Hnx) • VT(ni)] [a(n 2 ) • VT(n 2 )] T(n 3 )> 



(14) 



cancels with Eq. ( 13 ) and its cyclic permutations in the 



flat-sky limit if the three points are closely separated rel- 
ative to the correlation length of the lensing deflection. 
To calculate terms to higher order in Cf is likely best 
done via the real-space three-point function, extending 
the two-point methods of Ref. [32]. However, a more 
accurate calculation of the lensed bispectrum seems un- 
necessary given the already small effect of O(Cf^) shape 
changes by lensing on /jyx constraints (see below). 

We now turn to the effects of lensing on the estima- 
tion of Jnl- This involves an aggregate projection of 
the observed bispectrum onto our bispectrum of interest. 
In this case, the shape changes induced by lensing have 
the potential effectively to modify the normalization of 
the estimator. This may be investigated analytically if 
one has the capability to calculate a complete set of bis- 
pectrum coefficients, however the computational require- 
ments of this are daunting, scaling naively as C(^f nax ). 
Rather than perform a direct calculation, we place limits 



on any modification to the normalization using our simu- 
lations. We find that the estimator normalization in the 
presence of lensing is given by (/]^£ al ) = (1.0005)/]^ for 
the local shape and ( fff L ) — (0.965)/^ for the equilat- 
eral shape. Both of these small modifications will only 
become relevant if a high-significance detection of non- 
Gaussianity is made for which their effect becomes com- 
parable to the statistical errors. The smoothing of the 
bispectrum which lensing induces does not strongly af- 
fect estimators such as /jvl which average over £. We 
therefore conclude that the shape changes due to lensing 
have a negligible effect on the normalization of both cur- 
rent and near-future estimates of Jnl- Note, however, 
that the lensed CMB power spectrum should be used in 
the weights in the estimator [Eq. Q] to maintain near 
optimality [T3J [35J . As CMB experiments push to higher 
multipoles, we expect that the additional power which 
lensing generates on small scales will begin to have more 
noticeable effects on the $nl normalization (although 
other sources of non-Gaussianity will also grow rapidly). 
For Planck, however, the non-Gaussian effect of lensing 
on the estimator normalization is small enough that it 
may be ignored. 

As a further technical point, when estimating the indi- 
vidual bispectrum components B^^ a from the Monte- 
Carlo simulations in this section, we reduced the mea- 
surement error from the Monte-Carlo samples by using 
the following estimator: 
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where ap m denotes a Monte-Carlo simulation with f* L = 
1 and af m denotes the Gaussian part of the simulated 
temperature field (with f^ L = 0). The angle brackets in 
Eq. ( 15 ) denote a Monte-Carlo average. For more details 



on the non-Gaussian simulations, see the appendix. Sub- 
tracting the Gaussian part is analogous to the trick used 
for the ISW-lensing bias in ^ III A [Eq. (10)]: the second 
term on the right-hand side of Eq. (15) has zero mean 
but improves the variance of the estimator. Similarly, 
we estimate the change in normalization of the estimator 
/nl due to lensing by using the estimator 
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denotes a Monte-Carlo simulation with f* L 



where an 

1 and no ISW-lensing correlation (i.e. Cj^ = 
C. Lensing- induced variance 



0). 



For analytical calculations of prospective Jnl sensitiv- 
ity, the CMB is usually assumed to be observed on the 
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Var(/£ 



Fisher Sim. no ISW Sim. with ISW 



Local 17.0 
Equilateral 3240 



18.7 ± 1.9 
3710 ± 430 



19.7 ± 1.9 
3720 ± 450 



TABLE II: Variance of /jvl estimates in the absence of pri- 
mordial non-Gaussianity for the local and equilateral shapes. 



full sky and perfectly Gaussian. In this case direct calcu- 
lation shows that the variance of the estimator in Eq. Q 
is given by: 



Here, C(Cj*) denotes a term proportional to the cross 
spectrum Cj^, and is the ISW-lensing bias studied in 
§111 A| Because the ISW-lensing bias has a negligible de- 
pendence on cosmological parameters within their cur- 
rently allowed regions, for the purposes of an analysis 
of primordial non-Gaussianity one can simply subtract 
the bias without changing the definition of the estimator 
Inl, provided that the bias has been reliably computed. 
We find that the ISW-lensing bias is quite accurately ap- 
proximated by the lowest-order approximation in Eq. (|9| . 
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F(B X ,B X ) 



(Gaussian statistics) 



(17) 



) term in Eq. ( 18 ) is due to the change 
in the observed bispectrum shape studied in { III B In 
contrast to an earlier analysis [25] , we find that the differ- 
ence between the unlensed and lensed bispectra is small, 
on the order of 10% for i < 2000, in both the low-order 



(Note that we have assumed that each bispectrum B x is approximation of Eq. (|TT|) and in our fully lensed Simula- 
estimated independently; if multiple bispectra are being 
estimated jointly then the variance would be given by 
F^ 1 (B X ,B X ), where -F -1 is the matrix inverse.) 

However, because lensing generates non-Gaussianity 
even if the initial conditions are Gaussian, the true vari- 



ance may be larger than the right-hand side of Eq. (17) 
due to the additional connected three-, four- and six- 
point terms that are introduced. This issue has been 
studied for power spectrum estimation, where the excess 
estimator variance is small for £ max = 2500 in tempera- 
ture [HJ H2] > although it can be important when estimat- 
ing the i?-mode power spectrum in polarization [43H45] . 

Using our lensed simulations from the previous sections 
with the primordial f^L = 0, we may fit the distribution 
of estimated f^L values assuming that /jvl itself is ap- 
proximately Gaussian distributed. We do this both for 
simulations with the fiducial Cj^, as well as for simula- 
tions with Cj^ — 0, to investigate how the cosmic vari- 
ance of the ISW-lensing correlation affects the estimator 
variance. For the simulations with nonzero Cj^, we sub- 
tract the analytically calculated ISW-lensing bias from 
the Jnl estimates. Our findings are given in Table [ll| 
For both the local and equilateral configurations, there 
is marginal evidence for an increased /tvl variance due 
to the non-Gaussianity induced by lensing. This small 
increase in the variance is not enough to affect signifi- 
cantly current Fisher estimates. In a realistic analysis 
of non-Gaussianity with a sky-cut and inhomogeneous 
noise, errors on /jvl are usually estimated from Monte- 
Carlo simulations with /jvl = 0, and this increased vari- 
ance will be incorporated automatically provided that 
the simulations are correctly lensed. 



IV. DISCUSSION 

In the presence of lensing, the expectation value of the 
estimator f^ L can be written in the heuristic form 



tions. The shape modification has an even smaller effect 
on the effective normalization of both the local and equi- 
lateral /jvl estimators, and should not be important for 
Planck given its forecasted sensitivity and assuming that 
any detected Jnl will lie within current observational 
bounds. For example, were non-Gaussianity with the lo- 
cal shape to be detected at = 60 ±5, the bias due to 
neglecting the lens-induced shape change would be only 
A/j^£ 0.03. Similarly, a detection of the equilateral 
shape with f^ L — 250 ± 60 would have a bias A/jy^ ps 9. 
It is fortunate that lensing of the primordial bispectrum 
can be ignored since this spoils the separability of the bis- 
pectrum of the local and equilateral models that is key 
to fast fpfi, analyses. 

Finally, the variance of the /jvl estimates may be anal- 
ogously written 



Var(/*J 



F(B X ,B X ) 



0(Cf 



(19) 



The first term is the variance that would be obtained 
if the CMB were a Gaussian field. The second term in 



Eq. ( 19 1 represents excess variance due to non-Gaussian 



(18) 



statistics of the lensed CMB; our simulations suggest 
that it is small for a Planck-like experiment. Therefore, 
current Fisher-matrix forecasts for non-Gaussianity con- 
straints from Planck may still be considered accurate. 
Furthermore, any excess variance will be automatically 
incorporated into any future analysis which uses lensed 
simulations to determine the estimator variance in the 
presence of a sky-cut and inhomogeneous noise. 

Gravitational lensing subtly modifies the observed 
CMB. The bias which the ISW-lensing correlation in- 
troduces will need to be accurately subtracted in the lo- 
cal model. Beyond this, however, the non-Gaussianities 
which lensing creates represent small enough perturba- 
tions to the observed CMB that they may be safely ne- 
glected at Planck resolution. 
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Appendix: Non-Gaussian simulations 



In this work we generate non-Gaussian simulations us- 
ing the algorithm of [Hj. We simulate a Gaussian real- 
ization af m from the power spectrum Ci, and the non- 



Gaussian part of Eq. i6h is then generated by 



NG 



B 



x 



t i2 e 3 

m m 2 m 3 



h 2 m 2 



G* 

l t 3 m 3 



(20) 



Note that a£ m has zero mean for £ ^ 0. Since we do 
not simulate the monopole it is not necessary to subtract 
the mean explicitly. The above algorithm is a completely 
general method to create weakly non-Gaussian simula- 
tions with specified power spectrum Cg and bispectrum 
B^ 2 ^3, provided that contributions of order 0(fff L ) and 
higher can be neglected. More precisely, the two-, three-, 
and connected TV-point functions of the simulated field 
satisfy 



{ a l 1 rn 1 a ^ra 2 ) 
(f^fimi ^£2^2 ^3^3 ) 



\mi m 2 m 3 / 
0{f^ 2 ), (JV>4). (21) 



This simulation method is very general and computation- 
ally efficient, but one caveat is that terms of 0(f% L ) and 
higher are not explicitly controlled. For many purposes 
this caveat is irrelevant; to consider a concrete exam- 
ple, suppose we have an estimator (£/Af) of f NL whose 
overall normalization Af is unknown, and we are using 
the simulations to determine Af. (This approach to nor- 
malizing the estimator was used to analyze WMAP data 
in [ini EI].) Then we can compute Af as the following 
Monte-Carlo average: 



Af 




(22) 



In this form, it is clear that the terms containing two 
or more powers of /jvl are irrelevant, since they will not 
contribute to the derivative on the right-hand side. How- 
ever, there are cases where these terms are important; a 
concrete example would be the 0(ffj L ) contribution to 
the variance, Var(/jvL), of the /jvl estimator which was 
studied in [JHl HZ] • This term is proportional to the con- 
nected four-point function, and is not simulated reliably 
by the simulation algorithm in Eq. ( 20 ) . 



For the local shape, we found that the "subleading" 
Cf contribution to the power spectrum in Eq. (21 1 is 



spuriously large on large scales when using this simula- 
tion algorithm (Fig. [3]). It is frequently possible to work 



2 Note that we have written the simulation algorithm in a 
harmonic-space form which appears to have computational cost 
C(^max)i DU t f° r almost all bispectra of theoretical interest, there 
is a n eq uivalent position-space form with cost C(£^ lax ); see e.g. 
Eq. p4) below for the case of the local shape. 
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(1) Exact 

(2) General 

(3) Modified 
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FIG. 3: Power spectrum of the non-Gaussian component 
a em' l0C m t ne local model from three different simulation tech- 
niques: (1) uses the "exact" simulation algorithm of [48]; (2) 
uses the general algori thm of Eq. ( 24 1 ; and (3) uses the mod- 
ified algorithm of Eq. (271. 



around this problem if one is aware of it. For exam- 
ple, where simulations are used to determine the esti- 
mator normalization Af [Eq. ( |22[ )], the machinery from 
[H] allows the derivative with respect to f^L to be com- 
puted exactly, with zero contribution from higher powers 
°f Inl- However, our preferred solution is to make a 
simple change to the simulation algorithm for the spe- 
cial case of the local shape, which will give a reasonable 
power spectrum on large scales, as we now explain. 

First, note that in the local model, the reduced bispec- 
trum can be written for = 1 as 



7 local r* 

<^1^3 _ Z 



dr r 2 a tl (r)& a (r)/? 4 (r) + eye, (23) 



where the functions cxi(r) and fti{r) involve integrals over 
wavenumber k of the CMB transfer functions, spherical 
Bessel functions, je(kr), and the primordial power spec- 
trum; see e.g. [33] for details. Using this in Eq. (20) 
gives 



NG.loc 
1 lm 



= / dr 



fc(r) ^ d 2 nY; m (n)A(r 7 n)B(r,n) 



l - ai {r) y d 2 nY; m (n)B(r,nf 



(24) 



where we have defined 

A(r,n) = J2 a dr)af m Y em (n)/C e 

em 

B(r,fi) = J2^( r ) a fmYem(n)/Ci 



(25) 
(26) 
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If we modify the simulation algorithm by defining 



a 



,JVG,loc' 
im 



J drr 2 a e (r) (J d 2 nY e * m {n)B(r,n) 



(27) 

then the power spectrum on large scales is much smaller 
(see Fig. |3| , while it is easy to see that the expressions 
for the TV-point functions in Eq. (21 1 are unmodified. 



For the local shape, there is a different simulation al- 
gorithm available [35] which is "exact" , in the sense that 
the iV-point functions are reliably simulated (including 
subleading terms) for arbitrary N. This algorithm works 
by simulating the initial Newtonian potential $(r) inside 
a spherical shell which contains the surface of last scat- 
tering, and is thick enough to include all points inside 
the causal horizon at recombination. In Fig. [3j we also 
show the power spectrum of non-Gaussian simulations 
obtained using the exact simulation algorithm. It is seen 
that the general simulation algorithm [Eq. ( |20| )] produces 
a large-scale power spectrum which is much larger than 
the exact algorithm, but our modified algorithm for the 
local shape [Eq. (27)] produces a power spectrum which 



is smaller than the exact algorithm by an 0(1) factor. 
To understand intuitively why the modification pro- 



posed for the local shape in Eq. (27 1 eliminates the spu- 



riously large power spectrum on large scales which is 
generated by the general algorithm [Eq. (20)], we note 
that the modified algorithm has the following maximum- 
likelihood interpretation. Starting from a Gaussian CMB 
realization af m , suppose that we "guess" the potential 
$(r) throughout the Hubble volume, by finding the po- 
tential which maximizes the likelihood 



exp 



d 3 k |$(k)p 
(2tt) 3 P(k) 



(28) 



subject to the constraint that the potential $(r) projects 
to the observed CMB af . Suppose we then define a 



by projecting the potential $(r) 2 to an observed set of 
CMB multipoles. A short calculation then shows that the 
operation af m — > $(r) — > defined by this procedure 
is the same as the modified algorithm defined in Eq. ( p7| . 
Thus our modified algorithm for the local shape is closely 
related to the exact algorithm, but it produces a some- 
what smaller power spectrum because the deprojection 
operation ai m — > $(r) only generates power in one ra- 
dial mode for each value of £. In contrast, the general 



simulation algorithm in Eq. ( 20 1 does not appear to have 



any direct physical correspondence with the exact algo- 
rithm, and there is no guarantee that the power spectra 
are comparable. 

For the equilateral shape, we find that the 0(ff {L ) 
component of the power spectrum is reasonable on all 
scales and no modified version of the general simulation 
algorithm seems to be necessary. Formulating an exact 
simulation algorithm for the equilateral shape (i.e. an 
algorithm which precisely simulates A-point correlations 
for A > 4) has not been done; this would presumably 
require going back to the inflationary physics. 

In summary, there are several possible non-Gaussian 
simulation algorithms, with trade-offs as follows. For the 
equilateral shape, the general algorithm in Eq. (20) is the 



only known simulation procedure, and does not appear 
to contain any problems such as a spuriously large power 
spectrum. For the local shape, there is an algorithm [15] 
which is exact but computationally expensive (roughly 
100 CPU hours for the resolution requirements of this 
paper). This algorithm must be used for studies which 
require 0(f%f L ) and higher contributions to the iV-point 



functions to be simulated precisely (such as the 0(f% L ) 
contribution to Var(/)^£) studied in [IoT,l47| ). Otherwise, 
the algorithm proposed here for the local shape [Eq. ( 27 )] 



NG 
(m 



provides an alternative to the exact algorithm which is 
closely related but rather faster (roughly 20 CPU min- 
utes). 



